#include <QtCore/QCoreApplication>
#include <QPluginLoader>
#include <QDir>
#include <QFile>

#include "../plugin.h"
#include "../data/fileplugin.h"
#include "../data/reg3dvolumedata.h"
#include "../contree/contreeconstructor3d.h"
#include "../contree/augmentedctree.h"
#include "../contree/ctreecontourextractor.h"
#include "../math/timer.h"
#include "../mact/dctconstructor.h"
#include "../mact/dualgraph.h"
#include "../shape/sphere.h"
#include "../contour/mcconextractor.h"
#include "../sdf/geometryflow.h"
#include "../sdf/distancetransform.h"

//#include <boost/config.hpp>
#include <iostream>
#include <vector>
#include <algorithm>
#include <utility>
//#include <boost/graph/adjacency_list.hpp>
//#include <boost/graph/connected_components.hpp>


using namespace PBVP;

class MaximumPoint {
public:
	MaximumPoint(float _x = .0f, float _y = .0f, float _z = .0f, float _f = .0f) {
		x = _x;
		y = _y;
		z = _z;
		f = _f;
	}

	float x, y, z;
	float f;

	friend bool operator < (const MaximumPoint& p1, const MaximumPoint& p2) {
		return p1.f < p2.f;
	}

	//bool operator < (const MaximumPoint& p) const {
	//	return f < p.f;
	//}
};


int main(int argc, char *argv[])
{

	QCoreApplication a(argc, argv);
	//return a.exec();

	// load plugins
	QDir pluginsDir = QDir(a.applicationDirPath());
	std::cout << pluginsDir.absolutePath().toStdString() << std::endl;
	if (pluginsDir.dirName().toLower() == "debug"
		|| pluginsDir.dirName().toLower() == "release") {
			pluginsDir.cdUp();
	}
	// pluginsDir.cdUp();
	pluginsDir.cd("plugins");

	QStringList pluginEntries = pluginsDir.entryList(QDir::Files);

	FilePluginInterface* rawivPlugin;
	FilePluginInterface* rawPlugin;
	FilePluginInterface* rawncPlugin;
	FilePluginInterface* pdbPlugin;

	// check additional potential plugin path if necessary
	if(pluginEntries.size() == 0) {
		pluginsDir = QDir(a.applicationDirPath());
		pluginsDir.cd("plugins");
		pluginEntries = pluginsDir.entryList(QDir::Files);
	}

	std::cout << "Loading plugins in dir" << pluginsDir.absolutePath().toStdString() << std::endl;

	foreach (QString fileName, pluginEntries) {
		if (fileName.contains("rawiv"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin) {
				rawivPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (rawivPlugin) {
					rawivPlugin->initialize(NULL);
					std::cout << rawivPlugin->description().toStdString() << std::endl;
				}
			}
		} else if (fileName.contains("rawfile"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin) {
				rawPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (rawPlugin) {
					rawPlugin->initialize(NULL);
					std::cout << rawPlugin->description().toStdString() << std::endl;
				}
			}
		} else if (fileName.contains("rawncfile"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin)
			{
				rawncPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (rawncPlugin)
				{
					rawncPlugin->initialize(NULL);
					std::cout << rawncPlugin->description().toStdString() << std::endl;
				}
			}
		} else if (fileName.contains("pdbfile"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin)
			{
				pdbPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (pdbPlugin)
				{
					pdbPlugin->initialize(NULL);
					std::cout << pdbPlugin->description().toStdString() << std::endl;
				}
			}
		}
	}

	Timer timer;
	timer.getCurrentTime();
	std::cout << "processing file " << argv[1] << std::endl;
	Reg3DVolumeData *pReg3Data = (Reg3DVolumeData *)(rawivPlugin->loadFile(argv[1]));

	// look for maxima in the volume function
	// For pocket function, we want f > 0 and no other maximum in
	// a sphere of radius f
	std::vector<MaximumPoint> maxvec;
	std::vector<MaximumPoint>::iterator it;
	unsigned int dim[3];
	pReg3Data->getDimension(dim);
	for(int k = 0; k < dim[2]; k++) {
		for(int j = 0; j < dim[1]; j++) {
			for(int i = 0; i < dim[0]; i++) {
				float f = pReg3Data->getValue(i, j, k);
				bool maxflag = true;
				for(int kk = std::max<int>(0, k-1); kk < std::min<int>(k+1, dim[2]); kk++) {
					for(int jj = std::max<int>(0, j-1); jj < std::min<int>(j+1, dim[1]); jj++) {
						for(int ii = std::max<int>(0, i-1); ii < std::min<int>(i+1, dim[0]); ii++) {
							if(ii != i || jj != j || kk != k) {
								float f2 = pReg3Data->getValue(ii, jj, kk);
								if(f2 >= f) maxflag = false;
							}	// if statement
						}
					}
				}
				if(maxflag && f > 0.5) {	// consider maximum points with SDF > 0.5
					float coord[3];
					pReg3Data->getCoordinates(coord, i, j, k);
					MaximumPoint mp(coord[0], coord[1], coord[2], f);
					maxvec.push_back(mp);

				}
			} // for loops
		}
	}
	std::sort(maxvec.begin(), maxvec.end());
	//for(int i = 0; i < maxvec.size(); i++) {
	//	std::cout << "MP #" << i << ": " << maxvec[i].x << " " << maxvec[i].y << " " << maxvec[i].z << " " << maxvec[i].f << std::endl;
	//}
	int count = 0;
	GeometryData *pSpheres = new GeometryData();
	int size = maxvec.size();
	/*typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Graph;
	Graph G(size);
	for(int i = size-1; i >= 0; i--) {
		for(int j = i-1; j >= 0; j--) {
			MaximumPoint p1 = maxvec[i];
			Vector3 v1(p1.x, p1.y, p1.z);
			MaximumPoint p2 = maxvec[j];
			Vector3 v2(p2.x, p2.y, p2.z);
			if(v1.distance(v2) < (p1.f + p2.f)) {
				boost::add_edge(i, j, G);
			}
		}
	}*/
	/*std::vector<int> component(size);
	count = boost::connected_components(G, &component[0]);
	std::cout << "Total number of components: " << count << std::endl;
	int nc = component[size-1];		// component num of the large sphere
	for (int i = 0; i < size; i++)
	{
		if(component[i] == nc) { // spheres in the same component as the largest one
			MaximumPoint mp = maxvec[i];
			Vector3 v(mp.x, mp.y, mp.z);
			Sphere *s = Sphere::makeSphere(v, mp.f);
			GeometryData *pGeo = s->toGeometry(0);
			pSpheres->merge(pGeo);
			delete pGeo;
			delete s;
		}
	}*/

	// consider the pocket containing the largest sphere now
	MaximumPoint mp = maxvec.back();	// the largest sphere
	Vector3 v(mp.x, mp.y, mp.z);
	std::vector<MaximumPoint> pocket1;
	pocket1.push_back(mp);
	for (int i = 0; i < size; i++)
	{
		MaximumPoint p2 = maxvec[i];
		Vector3 v2(p2.x, p2.y, p2.z);
		if(v.distance(v2) < (mp.f + p2.f)) {
			pocket1.push_back(p2);
			//Sphere *s = Sphere::makeSphere(v2, p2.f);
			//GeometryData *pGeo = s->toGeometry(0);
			//pSpheres->merge(pGeo);
			//delete pGeo;
			//delete s;
		}
	}
	// construct the union of the spheres in the largest pocket
	float orig[3], span[3];
	Reg3DVolumeData *pSphVol;
	pReg3Data->getOrigin(orig);
	pReg3Data->getSpan(span);
	printf("span = (%f %f %f)\n", span[0], span[1], span[2]);

	/*pSphVol = new Reg3DVolumeData(dim);
	pSphVol->setOrigin(orig);
	pSphVol->setSpan(span);
	pSphVol->setNumOfVariables(1);
	float *sphdata = new float[dim[0]*dim[1]*dim[2]];
	for (int k = 0; k < dim[2]; k++)
	{
		for (int j = 0; j < dim[1]; j++)
		{
			for (int i = 0; i < dim[0]; i++)
			{
				int nsp = pocket1.size();
				Vector3 grid(orig[0]+span[0]*i, orig[1]+span[1]*j, orig[2]+span[2]*k);
				Vector3 cen(pocket1[0].x, pocket1[0].y, pocket1[0].z);
				int idx = k*dim[1]*dim[0] + j*dim[0] + i;
				sphdata[idx] = grid.distance(cen) - pocket1[0].f;
				for(int n = 1; n < nsp; n++) {
					Vector3 cp(pocket1[n].x, pocket1[n].y, pocket1[n].z);
					float dis = grid.distance(cp) - pocket1[n].f;
					if(dis < sphdata[idx]) sphdata[idx] = dis;
				}
			}
		}
	}
	pSphVol->setDataArray(sphdata);*/

	GeometryData *pPocket =  (GeometryData *)(rawPlugin->loadFile("1C3G.pdb_pk1.raw"));
	DistanceTransform dst(pPocket);
	dst.setDimension(dim);
	double minExt[3], maxExt[3];
	pReg3Data->getMinMaxExt(minExt, maxExt);
	dst.setMinMaxExt(minExt, maxExt);
	pSphVol = dst.computeSDF();
	for (int i = 0; i < dim[0]*dim[1]*dim[2]; i++)
	{
		pSphVol->setValue(-pSphVol->getValue(i), i);
	}
	// compute the zero contour of the volume function of spheres
	MCContourExtractor conExtractor;
	pSpheres = conExtractor.extractContour(pSphVol, 0.2);
	rawPlugin->saveFile(pSpheres, "spheres.raw");
	rawivPlugin->saveFile(pSphVol, "spheres.rawiv");
	delete pSpheres;
	// square the pocket SDF
	int nxyz = dim[0]*dim[1]*dim[2];
	for (int i = 0; i < nxyz; i++)
	{
		float val = pReg3Data->getValue(i)*pReg3Data->getValue(i);
		// float val = fabs(pReg3Data->getValue(i));
		pReg3Data->setValue(val, i);
	}
	GeometryFlow flow(pSphVol, pReg3Data);
	for (int i = 0; i < 50; i++)
	{
		flow.evolve(1);
		printf("Iteration # %d\n", i+1);

		pSpheres = conExtractor.extractContour(pSphVol, 0);
		char fname[256];
		sprintf(fname, "spheres%d.raw\0", i+1);
		rawPlugin->saveFile(pSpheres, fname);
		// sprintf(fname, "spheres%d.rawiv\0", i+1);
		// rawivPlugin->saveFile(pSphVol, fname);

		if (i > 0 && i%10 == 0)
		{
			flow.reinitilization();
		}
		delete pSpheres;
	}

	/*while(!maxvec.empty()) {
		MaximumPoint mp = maxvec.back();
		maxvec.pop_back();
		count ++;
		Vector3 v(mp.x, mp.y, mp.z);

		it = maxvec.begin();
		while(it != maxvec.end()) {
			MaximumPoint p2 = *(it);
			Vector3 v2(p2.x, p2.y, p2.z);
			if(v.distance(v2) < (mp.f + p2.f)) {
				it = maxvec.erase(it);
			} else {
				++it;
			}
		}
		if(count <= 100) {  // the first 100 non-overlapping spheres
			Sphere *s = Sphere::makeSphere(v, mp.f);
			GeometryData *pGeo = s->toGeometry(0);
			pSpheres->merge(pGeo);
			delete pGeo;
			delete s;
		}
		std::cout << "CP #" << count <<": " << mp.x << " " << mp.y << " " << mp.z << " " << mp.f << std::endl;
	}*/

	std::cout << "Elapsed time = " << timer.getElapsedTime() << std::endl;
	delete pReg3Data;
}
